load('simulation.RData')  


library(tidyverse)
## plot the sample size requirement
sample.size.DE.n20 <- data.frame(ICC = r.set, sample.size= J.DE.n20, n = "20", effect = "DE")
sample.size.DE.n100 <- data.frame(ICC = r.set, sample.size= J.DE.n100, n = "100", effect = "DE")
sample.size.MDE.n20 <- data.frame(ICC = r.set, sample.size= J.MDE.n20, n = "20", effect = "MDE")
sample.size.MDE.n100 <- data.frame(ICC = r.set, sample.size= J.MDE.n100, n = "100", effect = "MDE")
sample.size.SE.n20 <- data.frame(ICC = r.set, sample.size= J.SE.n20, n = "20", effect = "SE")
sample.size.SE.n100 <- data.frame(ICC = r.set, sample.size= J.SE.n100, n = "100", effect = "SE")

sample.size <- rbind(sample.size.DE.n20, 
                     sample.size.DE.n100,
                     sample.size.MDE.n20, 
                     sample.size.MDE.n100,
                     sample.size.SE.n20, 
                     sample.size.SE.n100)

label.effect <- c("DE" = "Direct Effect", "MDE" ="Marginal Direct Effect","SE"= "Spillover Effect")
ggplot(data = sample.size)+
  geom_line(mapping = aes(x = ICC, y =sample.size, linetype = n),show.legend = FALSE)+
  labs(x="Intracluster correlation coeffcient", y="Number of clusters")+
  scale_x_continuous(limits=c(0.1,1))+
  facet_wrap(~effect, labeller = labeller(effect = label.effect))+
  theme_bw()

ggsave("/Users/Zhichao/Dropbox/github/kosuke/mismatch/interference2/paper/figs/samplesize.pdf",height=4,width=8)

###  equal cluster size
p1.DE.n20.rho1 <- data.frame(ICC = r.set, power = power.DE.n20.equal.rho0, rho="0",n = "20", effect = " 'Direct Effect' ")
p1.DE.n20.rho2 <- data.frame(ICC = r.set, power = power.DE.n20.equal.rho03, rho="0.3",n = "20", effect = " 'Direct Effect' ")
p1.DE.n20.rho3 <- data.frame(ICC = r.set, power = power.DE.n20.equal.rho06, rho="0.6",n = "20", effect = " 'Direct Effect' ")
p1.MDE.n20.rho1 <- data.frame(ICC = r.set, power = power.MDE.n20.equal.rho0, rho="0",n = "20", effect = " 'Marginal Direct Effect' ")
p1.MDE.n20.rho2 <- data.frame(ICC = r.set, power = power.MDE.n20.equal.rho03, rho="0.3",n = "20", effect = " 'Marginal Direct Effect' ")
p1.MDE.n20.rho3 <- data.frame(ICC = r.set, power = power.MDE.n20.equal.rho06, rho="0.6",n = "20", effect = " 'Marginal Direct Effect' ")
p1.SE.n20.rho1 <- data.frame(ICC = r.set, power = power.SE.n20.equal.rho0, rho="0",n = "20", effect = " 'Spillover Effect' ")
p1.SE.n20.rho2 <- data.frame(ICC = r.set, power = power.SE.n20.equal.rho03, rho="0.3",n = "20", effect = " 'Spillover Effect' ")
p1.SE.n20.rho3 <- data.frame(ICC = r.set, power = power.SE.n20.equal.rho06, rho="0.6",n = "20", effect = " 'Spillover Effect' ")

p1.DE.n100.rho1 <- data.frame(ICC = r.set, power = power.DE.n100.equal.rho0, rho="0",n = "100", effect = " 'Direct Effect' ")
p1.DE.n100.rho2 <- data.frame(ICC = r.set, power = power.DE.n100.equal.rho03, rho="0.3",n = "100", effect = " 'Direct Effect' ")
p1.DE.n100.rho3 <- data.frame(ICC = r.set, power = power.DE.n100.equal.rho06, rho="0.6",n = "100", effect = " 'Direct Effect' ")
p1.MDE.n100.rho1 <- data.frame(ICC = r.set, power = power.MDE.n100.equal.rho0, rho="0",n = "100", effect = " 'Marginal Direct Effect' ")
p1.MDE.n100.rho2 <- data.frame(ICC = r.set, power = power.MDE.n100.equal.rho03, rho="0.3",n = "100", effect = " 'Marginal Direct Effect' ")
p1.MDE.n100.rho3 <- data.frame(ICC = r.set, power = power.MDE.n100.equal.rho06, rho="0.6",n = "100", effect = " 'Marginal Direct Effect' ")
p1.SE.n100.rho1 <- data.frame(ICC = r.set, power = power.SE.n100.equal.rho0, rho="0",n = "100", effect = " 'Spillover Effect' ")
p1.SE.n100.rho2 <- data.frame(ICC = r.set, power = power.SE.n100.equal.rho03, rho="0.3",n = "100", effect = " 'Spillover Effect' ")
p1.SE.n100.rho3 <- data.frame(ICC = r.set, power = power.SE.n100.equal.rho06, rho="0.6",n = "100", effect = " 'Spillover Effect' ")


power1 <- rbind(p1.DE.n20.rho1,
                p1.DE.n20.rho2,
                p1.DE.n20.rho3,
                p1.MDE.n20.rho1,
                p1.MDE.n20.rho2,
                p1.MDE.n20.rho3,
                p1.SE.n20.rho1,
                p1.SE.n20.rho2,
                p1.SE.n20.rho3,
                p1.DE.n100.rho1,
                p1.DE.n100.rho2,
                p1.DE.n100.rho3,
                p1.MDE.n100.rho1,
                p1.MDE.n100.rho2,
                p1.MDE.n100.rho3,
                p1.SE.n100.rho1,
                p1.SE.n100.rho2,
                p1.SE.n100.rho3
                )

power1$rho <- factor(power1$rho, levels = c("0","0.3","0.6"),
                  labels=c(expression(paste(rho, "=0")),expression(paste(rho, "=0.3")),expression(paste(rho, "=0.6"))  )
)
ggplot(data = power1)+
  geom_line(mapping = aes(x = ICC, y =power, linetype = n),show.legend = FALSE)+
  geom_hline(yintercept = 0.8,linetype = 2,color="grey")+
  labs(x="Intracluster correlation coeffcient", y="Power")+
  facet_grid(rho~effect,labeller=label_parsed)+
  scale_y_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1),limits=c(0,1))+
  scale_x_continuous(limits=c(0.1,1))+
  theme_bw()

ggsave("power1.pdf",height=4,width=8)




#### unequal cluster size

p2.DE.n20.rho1 <- data.frame(ICC = r.set, power = power.DE.n20.unequal.rho0, rho="0",n = "20", effect = " 'Direct Effect' ")
p2.DE.n20.rho2 <- data.frame(ICC = r.set, power = power.DE.n20.unequal.rho03, rho="0.3",n = "20", effect = " 'Direct Effect' ")
p2.DE.n20.rho3 <- data.frame(ICC = r.set, power = power.DE.n20.unequal.rho06, rho="0.6",n = "20", effect = " 'Direct Effect' ")
p2.MDE.n20.rho1 <- data.frame(ICC = r.set, power = power.MDE.n20.unequal.rho0, rho="0",n = "20", effect = " 'Marginal Direct Effect' ")
p2.MDE.n20.rho2 <- data.frame(ICC = r.set, power = power.MDE.n20.unequal.rho03, rho="0.3",n = "20", effect = " 'Marginal Direct Effect' ")
p2.MDE.n20.rho3 <- data.frame(ICC = r.set, power = power.MDE.n20.unequal.rho06, rho="0.6",n = "20", effect = " 'Marginal Direct Effect' ")
p2.SE.n20.rho1 <- data.frame(ICC = r.set, power = power.SE.n20.unequal.rho0, rho="0",n = "20", effect = " 'Spillover Effect' ")
p2.SE.n20.rho2 <- data.frame(ICC = r.set, power = power.SE.n20.unequal.rho03, rho="0.3",n = "20", effect = " 'Spillover Effect' ")
p2.SE.n20.rho3 <- data.frame(ICC = r.set, power = power.SE.n20.unequal.rho06, rho="0.6",n = "20", effect = " 'Spillover Effect' ")

p2.DE.n100.rho1 <- data.frame(ICC = r.set, power = power.DE.n100.unequal.rho0, rho="0",n = "100", effect = " 'Direct Effect' ")
p2.DE.n100.rho2 <- data.frame(ICC = r.set, power = power.DE.n100.unequal.rho03, rho="0.3",n = "100", effect = " 'Direct Effect' ")
p2.DE.n100.rho3 <- data.frame(ICC = r.set, power = power.DE.n100.unequal.rho06, rho="0.6",n = "100", effect = " 'Direct Effect' ")
p2.MDE.n100.rho1 <- data.frame(ICC = r.set, power = power.MDE.n100.unequal.rho0, rho="0",n = "100", effect = " 'Marginal Direct Effect' ")
p2.MDE.n100.rho2 <- data.frame(ICC = r.set, power = power.MDE.n100.unequal.rho03, rho="0.3",n = "100", effect = " 'Marginal Direct Effect' ")
p2.MDE.n100.rho3 <- data.frame(ICC = r.set, power = power.MDE.n100.unequal.rho06, rho="0.6",n = "100", effect = " 'Marginal Direct Effect' ")
p2.SE.n100.rho1 <- data.frame(ICC = r.set, power = power.SE.n100.unequal.rho0, rho="0",n = "100", effect = " 'Spillover Effect' ")
p2.SE.n100.rho2 <- data.frame(ICC = r.set, power = power.SE.n100.unequal.rho03, rho="0.3",n = "100", effect = " 'Spillover Effect' ")
p2.SE.n100.rho3 <- data.frame(ICC = r.set, power = power.SE.n100.unequal.rho06, rho="0.6",n = "100", effect = " 'Spillover Effect' ")



power2 <- rbind(p2.DE.n20.rho1,
                p2.DE.n20.rho2,
                p2.DE.n20.rho3,
                p2.MDE.n20.rho1,
                p2.MDE.n20.rho2,
                p2.MDE.n20.rho3,
                p2.SE.n20.rho1,
                p2.SE.n20.rho2,
                p2.SE.n20.rho3,
                p2.DE.n100.rho1,
                p2.DE.n100.rho2,
                p2.DE.n100.rho3,
                p2.MDE.n100.rho1,
                p2.MDE.n100.rho2,
                p2.MDE.n100.rho3,
                p2.SE.n100.rho1,
                p2.SE.n100.rho2,
                p2.SE.n100.rho3
)


power2$rho <- factor(power2$rho, levels = c("0","0.3","0.6"),
                     labels=c(expression(paste(rho, "=0")),expression(paste(rho, "=0.3")),expression(paste(rho, "=0.6"))  )
)

ggplot(data = power2)+
  geom_line(mapping = aes(x = ICC, y = power, linetype = n),show.legend = FALSE)+
  geom_hline(yintercept = 0.8,linetype = 2,color="grey")+
  labs(x="Intracluster correlation coeffcient", y="Power")+
  facet_grid(rho~effect,labeller=label_parsed)+
  scale_y_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1),limits=c(0,1))+
  scale_x_continuous(limits=c(0.1,1))+
  theme_bw()

ggsave("power2.pdf",height=4,width=8)


